function [parGMM_2, fGMM_2, retcodeGMM_2,parGMM_1, fGMM_1, retcodeGMM_1] = GMMmultiloss3(par0);
global P ERROR_SERIES INSTRUMENTS N_ITER WEIGHT_WT_1 WEIGHT_WT_2 PAR_N_ITER_1 PAR_N_ITER_2;

% GMM_MULTI_LOSS_2.M            Version I.K. 22/12/2006
%   Constructs the GMM estimates for multivariate loss parameters in a bivariate loss:
%       L2(p,alpha,e) = 0.5*(norm(e,p)+t'*e)*norm(e,p)^(p-1)
% and a sum of univatiate losses: 
%       L1(p,alpha_1,e_1) + ... + L1(p,alpha_d,e_d) = (alpha - t.*(e < 0))'*(abs(e).^p)
% where
%   forecast error e = (e_1,...,e_d)'
%   exponent P (GLOBAL variable) is scalar AND FIXED!!!
%   asymmetry vector t = (2*alpha_1-1,...,2*alpha_d-1)'
%  
% Inputs:   ERROR_SERIES = time series of forecast error N-vectors (GLOBAL variable)
%           INSTRUMENTS = time series of K-vector instruments (GLOBAL variable)
%           par0 = (N+1)-vector of starting parameter values, par0 = (p0; alpha0)
% Output:   parGMM = (max_iter x (N+1)) matrix of GMM estimates, parGMM(n,:) = (p; alpha')'
%           fGMM = (max_iter x 1) vector of optimal values of the GMM objective function 
%           retcodeGMM = (max_iter x 1) vector of CV codes
% 
% COMMENTS:
%   We minimize the objective function Q(par,Yt) and use the notations in Hamilton, ch 14.

options = optimset('MaxIter',200,'MaxFunEvals',1000,'TolFun', 1e-8,'Display','off',...
    'LargeScale','off','GradObj','off');

% Linear inequality constraints A.x <= b on the parameters of the model
delta = 0.0001; % some small number
[T,N] = size(ERROR_SERIES);
A = zeros(2*N,N);
b = [];
for i=1:1:N
    A(2*i-1,i) = 1;       % alpha_i < 1
    A(2*i,i) = -1;        % alpha_i > 0
    b = [b;[1-delta; -delta]];
end

% Recursive estimation 
max_iter = 50;
q = floor(1+T^0.2);     % cf Hamilton eq 14.1.19
parGMM_1 = []; fGMM_1 = []; retcodeGMM_1 = [];
parGMM_2 = []; fGMM_2 = []; retcodeGMM_2 = [];
N_ITER = 1;
PAR_N_ITER_1 = par0;
PAR_N_ITER_2 = par0;
while N_ITER <= max_iter
    % Display N_ITER
    N_ITER;
    if N_ITER > 1
        WEIGHT_WT_1 = WeightWT_1(q,PAR_N_ITER_1);
        WEIGHT_WT_2 = WeightWT_2(q,PAR_N_ITER_2);
    end
    % multivariate loss
    [parGMM_N_ITER,fGMM_N_ITER,retcodeGMM_N_ITER,output_N_ITER,lambda_N_ITER,grGMM_N_ITER,hessGMM_N_ITER] = ...
        fmincon(@ObjFctGMMmulti2,PAR_N_ITER_2,A,b,[],[],[],[],[],options);
    PAR_N_ITER_2 = parGMM_N_ITER;
    parGMM_2 = [parGMM_2; parGMM_N_ITER'];
    fGMM_2 = [fGMM_2; fGMM_N_ITER];
    retcodeGMM_2 = [retcodeGMM_2; retcodeGMM_N_ITER];
    % naive sum of univariate losses
    [parGMM_N_ITER,fGMM_N_ITER,retcodeGMM_N_ITER,output_N_ITER,lambda_N_ITER,grGMM_N_ITER,hessGMM_N_ITER] = ...
        fmincon(@ObjFctGMMuni1,PAR_N_ITER_1,A,b,[],[],[],[],[],options);
    PAR_N_ITER_1 = parGMM_N_ITER;
    parGMM_1 = [parGMM_1; parGMM_N_ITER'];
    fGMM_1 = [fGMM_1; fGMM_N_ITER];
    retcodeGMM_1 = [retcodeGMM_1; retcodeGMM_N_ITER];
    N_ITER = N_ITER+1;
end


